home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / cfft.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  2KB  |  113 lines
  1. /*
  2.  *    Definitions for the double precision Complex FFT server
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)cfft.cc    2.1    8/18/89
  23.  */
  24.  
  25.  
  26.  
  27. #include "rw/DComplexFFT.h"
  28. #include "rw/mathpack.h"
  29.  
  30. DComplexFFTServer::DComplexFFTServer()
  31. {
  32.   server_N = 0;
  33. }
  34.  
  35. DComplexFFTServer::DComplexFFTServer(unsigned N)
  36. {
  37.   setOrder(N);
  38. }
  39.  
  40. DComplexFFTServer::DComplexFFTServer(const DComplexFFTServer& s)
  41.      : the_weights(s.the_weights)
  42. {
  43.   server_N = s.server_N;
  44. }
  45.  
  46. void
  47. DComplexFFTServer::operator=(const DComplexFFTServer& s)
  48. {
  49.   server_N = s.server_N;
  50.   the_weights.reference(s.the_weights);
  51. }
  52.  
  53. void
  54. DComplexFFTServer::setOrder(unsigned N)
  55. {
  56.   the_weights.resize(4*N+15);
  57.   fortran_int local = server_N = N;    // Local copy, converted to int
  58.   DCffti_(&local, the_weights.data());    // Call the NCAR weight routine
  59. }
  60.  
  61. /*
  62. Use the NCAR FFT routines.  See notes in class description header. 
  63. */
  64.  
  65. DComplexVec
  66. DComplexFFTServer::fourier(const DComplexVec& v)
  67. {
  68.   unsigned N = v.length();
  69.   if( N != server_N )setOrder(N);
  70.  
  71.   fortran_int temp_length = N;
  72.  
  73. #if COMPLEX_PACKS
  74.   /* Make a copy to be transformed: */
  75.   DComplexVec inverse = v.copy();
  76.   DCfftf_(&temp_length, inverse.data(), the_weights.data());
  77. #else
  78. You lose.    Code requires that sizeof(complex)==2*sizeof(double)
  79. #endif
  80.  
  81.   return inverse;
  82. }
  83.  
  84. /* Inverse transform */
  85.  
  86. DComplexVec
  87. DComplexFFTServer::ifourier(const DComplexVec& v)
  88. {
  89.   unsigned N = v.length();
  90.   if( N != server_N )setOrder(N);
  91.  
  92.   fortran_int temp_length = N;
  93.  
  94. #if COMPLEX_PACKS
  95.   /* Make a copy to be transformed: */
  96.   DComplexVec inverse = v.copy();
  97.   DCfftb_(&temp_length, inverse.data(), the_weights.data());
  98. #else
  99. You lose.    Code requires that sizeof(complex)==2*sizeof(double)
  100. #endif
  101.  
  102.   return inverse;
  103. }
  104.  
  105. /************  Useful utilities **************/
  106.  
  107. double
  108. spectralVariance(const DComplexVec& v)
  109. {
  110.   // Do not count the mean v(0):
  111.   return sum(norm(v.slice(1,v.length()-1,1)));
  112. }
  113.